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Type la Supernovae: An Asymmetric Deflagration Model 
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ABSTRACT 

We present the first high-resolution three-dimensional simulations of the de- 
flagration phase of Type la supernovae that treat the entire massive white dwarf. 
We report the results of simulations in which ignition of the nuclear burning oc- 
curs slightly off-center. The subsequent evolution of the nuclear burning is sur- 
prisingly asymmetric with a growing bubble of hot ash rapidly rising to the stellar 
surface. Upon reaching the surface, the mass of burned material is ~ O.O75M 
and the kinetic energy is 4.3 x 10 49 ergs. The velocity of the top of the rising 
bubble approaches 8000 kms" 1 . The amount of the asymmetry found in the 
model offers a natural explanation for the observed diversity in Type la super- 
novae. Our study strongly disfavors the classic central-ignition pure deflagration 
scenario by showing that the result is highly sensitive to details of the initial 
conditions. 

Subject headings: hydrodynamics — instabilities - stars:interior — supernovae:general 
— white dwarfs 



1. Introduction 

Type la supernovae are one class of stellar explosions noted for their extreme brightness 
and distinguished by the lack of hydrogen in their spectra. These events are thought to 
originate in close binary systems in which the primary component is a degenerate C-0 rich 
white dwarf (Hoyle & Fowler 1960) gaining mass accreted from a low- mass main sequence 
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companion. A thermonuclear explosion is believed to be initiated in the core of the white 
dwarf when its mass approaches the Chandrasekhar limit (Whelan & Iben 1973; Livio 2001). 
The exact nature of the ignition and subsequent evolution, however, is the subject of ongoing 
research. 

Regardless of the exact nature of the explosion, the nuclear energy released by the decay 
of 56 Ni and 56 Co makes Type la supernovae extremely bright optical events. In addition, 
observations indicate that Type la supernova light curves form a homogeneous class with 
relatively small intrinsic scatter (Filippenko 1997). The brightness and homogeneity of the 
light curves have earned these supernovae "standard candle" status for determining distances 
in the Universe and thus cosmological parameters (Sandage & Tammann 1993; Branch 1998; 
Garnavich et al. 1998; Perlmutter et al. 1999; Turner 2001; Tonry et al. 2003; Knop et al. 
2003, and references therein). Despite the apparent homogeneity of these events, some 
observations show exceptions (Li et al. 2001). 

In the classic Type la scenario, a carbon flash evolves into a large-scale explosion, 
eventually consuming the entire white dwarf. The thermonuclear flame propagates through 
the C-0 fuel as either a subsonic deflagration front (Nomoto, Sugimoto, & Neo 1976; Nomoto, 
Thielemann, & Yokoi 1984; Reinecke, Hillebrandt, & Niemeyer 2002; Gamezo et al. 2003) 
or a supersonic detonation wave (Arnett 1969; Boisseau et al. 1996) and releases sufficient 
energy to unbind the star. Models involving either a pure deflagration or a pure detonation, 
however, are unable to provide an explanation for both the observed expansion velocities and 
the spectrum produced by ejecta that are rich in intermediate-mass and iron-peak elements 
(Truran & Cameron 1971; Woosley & Weaver 1986; Woosley 1990). 

Considerable effort has gone into resolving the discrepancy with observations. One 
possibility is a transition of the initial deflagration front into a detonation as the burning 
front propagates outward through the outer layers of the white dwarf (Blinnikov & Khokhlov 
1987; Khokhlov 1991b; Niemeyer & Woosley 1997). Another possibility is that the initial 
deflagration dies out as a result of the expansion of the outer layers of the white dwarf. When 
these gravitationally-bound layers recollapse onto the white dwarf, a detonation ensues, 
incinerating the remaining nuclear fuel (Khokhlov 1991a). These models are often referred 
to as "delayed detonation" models. In addition, our knowledge regarding initial conditions 
characterizing the stellar core just prior to ignition is severely limited. Analytic models 
suggest rather complex flow patterns (Woosley, Wunsch & Kuhlen 2003), while the first 
multidimensional hydrodynamic models of stellar interiors became available only recently 
(Hdflich & Stein 2002). 

Despite these deficiencies, delayed detonation models of massive C-0 rich white dwarfs 
are capable of accounting for most Type la characteristics, including the observed expansion 
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velocities of silicon-group and iron-group nuclei. In this Letter we report on research into 
the initial deflagration phase of the explosion and present the results of simulations in which 
ignition occurs slightly off-center. 

2. Numerical model 

Modeling complex events like Type la supernovae is limited by available computing 
resources due to the vast range of length scales involved. Theoretical models of thermonuclear 
flame fronts indicate that, for conditions typical of white dwarf interiors, the flame thickness 
is between 8 and 12 orders of magnitude smaller than the stellar radius ( Timmes & Woosley 
1992). Because of this difference in length scales, multidimensional Type la models must 
make use of an appropriate subgrid model for the evolution of the thermonuclear burning 
front. 

A subsonic burning front originating near the center of a massive white dwarf is subject 
to several fluid instabilities (Niemeyer & Kerstein 1997). These instabilities, dominated 
by the Rayleigh- Taylor instability (RTI) on the largest scales, increase the surface area of 
the flame and its effective propagation speed. As was proposed by Khokhlov (1995), this 
increase in speed is offset by the effects of the RTI: the more quickly the flame advances, 
the more rapidly the bubbles produced by RTI merge. The surface area and flame speed, 
therefore, do not increase. Khokhlov's self-regulating mechanism was observed in numerical 
simulations of compressible as well as incompressible flows (Khokhlov 1995; Vladimirova 
& Rosner 2003), and current state-of-the-art numerical models (Reinecke, Hillebrandt, & 
Niemeyer 2002; Gamezo et al. 2003) indicate that this mechanism might be robust. We will 
report in an upcoming publication a series of simulations that confirm these results. 

In our modeling effort, we performed entire-star simulations of the deflagration phase 
of Type la supernovae using the FLASH code, a parallel adaptive-mesh multi-physics as- 
trophysical hydrodynamics code (Fryxell et al. 2000; Calder et al. 2002). The FLASH code 
solves the Euler equations for compressible flow and the Poisson equation for self-gravity. A 
recent addition is a custom implementation of the flame capturing scheme of Khokhlov (2001) 
in which the flame advances by evolving a passive scalar variable with an advection-reaction- 
diffusion equation. Reaction and diffusion coefficients are chosen to keep the interface several 
grid points wide and to propagate with the given speed. Specifically, we take the flame speed 
to be a maximum of a laminar flame speed and (Timmes & Woosley 1992) and a turbulent 
flame speed, based on the assumption that the turbulent burning on a macroscopic scales is 
driven by RTI Khokhlov (1995) and described above. The flame model controls the transi- 
tion from carbon to magnesium and produces the corresponding amount of energy. Further 
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transition to silicon and nickel are modelled as two consecutive stages of distributed burning 
with exponential decays on timescales that depend on temperature and density. 

For our initial configuration, we adopt a "cold" isothermal model of a white dwarf with 
a temperature of T wd = 5 x 10 7 K, a mass M wd = 1.36M , and a radius i? w( j = 2.13 x 10 8 cm. 
We found that this model was not in hydrostatic equilibrium after being interpolated onto 
the FLASH simulation mesh due to a small mismatch in discretization of the solution and 
slight differences in the equation of state between that of the initial analytic model and 
that used in the simulation. To bring the model into hydrostatic equilibrium, and therefore 
ensure that the subsequent evolution is not affected by this discrepancy, we used a variant of 
the relaxation method proposed by Arnett (1994). The remaining mismatch in the relaxed 
model produced velocities smaller than 10 kms -1 during test simulations executed for several 
seconds of evolution time. We will present these and other verification tests of the method 
in a forthcoming publication (Calder et al. 2004). 

For the multidimensional simulations, we mapped the relaxed model onto the adaptive 
simulation mesh generated according to the following refinement criteria. We forced the 
resolution in the innermost part of the grid containing the bulk of the star, r uni = 2300 
km, to be no worse than a predefined limit, Ar uni = 12.8 km. During the simulation, the 
same resolution was also enforced in any grid region where the density exceeded p un j = 
3 x 10 6 gem -3 . Grid refinement was also active in the region where local density variations 
were greater than 20%. We used the same criteria to resolve velocity variations when the 
velocity exceeds 100 kms" 1 . Finally, we resolved the region occupied by the flame front 
identified by the passive scalar at the highest resolution. These refinement criteria result in 
an initial grid configuration that uniformly resolves the interior of the star, and the resolution 
gradually decreases with radius outside of the star. During the evolution, the mesh adapts 
so that only the regions at high densities are refined. 

We initialized the flame front by placing a small spherical region of completely burned 
material in hydrostatic equilibrium with its surroundings near the stellar center. We note 
that our configuration differs significantly from an off-center ignition model considered 
by Niemeyer, Hillebrandt & Woosley (1996) in which the ignition region was placed much 
farther from the stellar center. We determined the subsequent evolution of the flame by 
integrating the hydrodynamic equations supplemented by the advection-diffusion-reaction 
equation. We used a tabular equation of state suitable for the dense stellar material and 
performed the simulation in Cartesian geometry in three dimensions. At the boundaries of 
the domain, we imposed hydrostatic outflow-only boundary conditions. We calculated the 
gravitational field assuming spherical symmetry, an assumption justified by the relatively 
small displacement of mass during the simulation. 
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3. Results 

In the first simulation, we used a computational domain with sides L = 6.5536 x 10 8 cm 
and an effective resolution of 1024 3 (6.4 km); in the second simulation, we used the same 
domain but with an effective resolution of 4096 3 (1.6 km). Due to the large size of the 
high-resolution simulation, the evolution after t = 0.76 s was calculated at an effective 
resolution of 3.2 km. We initiated the nuclear flame in a spherical region of radius of 50 km 
centered at (x, y,z) = (10,6,4) km (i.e., 12 km off-center). This displacement corresponds 
to about 0.25 of the radius of the ignition region, a few percent of the expected radius of 
the convective region in the center of the progenitor, and 6 x 10~ 3 of the radius of the star. 
This is significantly smaller than the 200 km displacement used by Niemeyer, Hillebrandt 
& Woosley (1996). Asymmetric placement of the ignition region was aimed at minimizing 
possible mesh discretization effects on the early flame evolution. (Verification tests with a 
centrally ignited flame indicated that mesh effects are strong enough to promote preferential 
instability growth along the grid axes.) 

We find that the evolution of the flame front (the bubble of hot ash) consists of three 
distinct phases. The first lasts for approximately 0.3-0.4 s and is a period of slow, mostly 
radial growth of the bubble. The position of the hot bubble changes only slightly, but one 
can observe a slow divergence of its shape from spherical symmetry, caused by the non- 
uniformity of the gravitational field. The next phase is characterized by substantial growth 
of the asymmetry created during the first phase and the development of RTI at the bubble 
surface. At this point, the difference between the turbulent flame speed at the '"top" (further 
away from the center of the star) and the "bottom" (closer to the center) of the bubble leads 
to further distortion of the bubble. Later in this phase, the top of the bubble becomes highly 
susceptible to RTI, which leads to the successive formation of a series of effusive plumules 
(see Fig. 1(a)). This structure rapidly grows in size and complexity (see Fig. 1(b) and Fig. 2), 
increasing the surface area of the flame and therefore the burning rate and rate of energy 
release (see Fig. 3). The whole structure starts to rapidly ascend toward the surface. As in 
the first phase, the bubble does not encounter significant density gradients, and the evolution 
of the surface is dominated by the flame evolution. 
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Fig. 1. — Flame surface evolution, a) The bubble at the moment the Rayleigh- Taylor 
unstable top surface begins to grow rapidly (i.e, at the end of the first phase of its evolution), 
b) The bubble as it approaches the stellar surface (i.e., at the end of the second phase 
of its evolution). The complex structure results from the merging of secondary and later 
generations of bubbles produced by RTI. Filaments and sheets of largely unburned material 
are present throughout the bubble. For comparison, the red sphere has a radius of 25 km in 
both images. 



Fig. 2. — Volume rendering of the bubble of hot ash as it approaches the stellar surface (i.e., 
at the end of the second phase of its evolution). The complex structure is formed by the 
merging of secondary and later generations of bubbles produced by RTI. The length scale in 
the image is identical to that of Fig. 1(b). 




Fig. 3. — Some global properties of the model as a function of time. Shown are the total 
binding energy (thick solid), the energy generation rate Qb (dotted), the flame front area 
A\y (thin), the total kinetic energy E^ (thick dashed), and the mass of the bubble (thick 
dotted). 
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The third phase of the evolution is marked by the time at which the bubble reaches 
the surface layers of the star where there is a large decreasing density gradient, and lasts 
for only ps 0.05 s. Here pressure effects become important to the evolution, and lateral 
expansion accompanies the forward motion. The evolution is extremely fast. At 0.78 s, the 
tip of the bubble is located at r ~ 1.51 x 10 s cm and by t = 0.83 s the bubble reaches a 
radius of ~ 1.88 x 10 8 cm. By the end of the simulation, the bubble has a radial extent 
of about 1500 km and a similar diameter. The interior of the bubble has a highly complex 
structure, with a tangled network of burned and partially burned material having bulk 
velocities approaching 10,000 kms -1 and lateral expansion in the range of 3000 km s" 1 . The 
rapid expansion of the bubble and the associated decrease of density and temperature result 
in a slowing of the thermonuclear reaction rates and flame extinction. The physics expected 
in the subsequent evolution are not properly accounted for in our model (any burning will 
no longer have the form of a well-defined deflagration front), so at this point we ended the 
simulations. 

Figure 3 depicts temporal evolution of select global quantities. The result is that ~ 
0.026M Q (« 2% of the stellar mass) burns into iron-peak elements. The rapidly rising 
bubble expels burned material originating in the center of the star with a characteristic 
velocity of 5000 kms -1 and with local velocities in excess of 10,000 kms -1 . Only about 
0.2% of the stellar mass has sufficient energy to leave the system, though. The rest of the 
burned material remains gravitationally bound and will fall back onto the surface of the star, 
thereby increasing the abundance of Ni-group elements in the surface layers. 



4. Discussion and Conclusions 

In this paper, we studied the hydrodynamics of thermonuclear deflagration of a massive 
white dwarf. We considered a flame ignited slightly off-center and followed its evolution in 
a full three-dimensional simulation of an entire white dwarf with no imposed symmetry. We 
observed the formation of a Rayleigh- Taylor-unstable bubble of hot nuclear ash and its rapid 
rise toward the stellar surface. The rising bubble takes the form of a typical mushroom- 
shaped plume and is composed of intermediate and iron-peak elements. As it rises, the 
bubble's motion is continuously powered by buoyancy due to heating from the energy released 
by nuclear burning, and supersonic speeds are achieved. Further, the morphology of the 
surface of the bubble, determined by RTI, displays characteristics similar to those obtained 
in earlier numerical studies of deflagrating white dwarfs. In particular, generations of merging 
bubbles form on the surface of the initial bubble. 

Although this model does not produce an immediate explosion, it provides a robust 
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mechanism for the transportation of heavy elements to the stellar surface. If combined 
with a successful explosion, such a surface composition might offer an explanation for the 
presence of iron-group elements at high velocities, as observed in SN 1991T-like supernovae 
(Filippenko 1997; Jeffrey et al. 1992). Our model produces ~ O.O3M of heavy elements 
at the stellar surface without invoking an ad hoc transition to detonation (Yamaoka et al. 
1992). 

Since the buoyancy- driven velocity created soon after the ignition quickly becomes com- 
parable to the expected velocity of the convective motions (Hofiich & Stein 2002), convection 
is unlikely to influence the bubble's evolution in any significant way. This conclusion is con- 
sistent with the results of Niemeyer, Hillebrandt & Woosley (1996). 

The fact that our initial conditions deviate from perfect symmetry by only a minute 
amount, but lead to a result drastically different from a central ignition model, indicates that 
central ignition models in general are unlikely to be typical initial conditions for the explosion. 
Our results therefore strongly disfavor the classic central-ignition pure deflagration scenario. 

In future studies, we plan to focus on the late time evolution of this model (Plewa 
et al. 2004). Recollapse and/or detonation are likely outcomes. Detailed calculations of 
the nucleosynthetic yields and further convergence studies will be required to provide a 
quantitative link between this model and observations. Finally, it will be important to study 
the evolution of the deflagration front in the presence of strong shear, which can influence 
the evolution of the rising bubble, in more detail. 

The authors thank Peter Hofiich and Alexei Khokhlov for their comments. The authors 
also acknowledge visualization support from the ANL Futures Laboratory. This work is 
supported in part by the U.S. Department of Energy under Grant No. B341495 to the 
Center for Astrophysical Thermonuclear Flashes at the University of Chicago. 
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